## Created by EmpowerStats @ Fri, 06 Dec 24 20:11:09 +0800## 
#******************** Regarding ALL Following R Code   ******************** 
#*** COPYRIGHT (c) 2010, 2021 X&Y Solutions, ALL RIGHT RESERVED *********** 
#******************* www.EmpowerStats.com ********************************* 
#************************************************************************** 
Sys.setlocale(category = 'LC_ALL', locale = 'English_United States.1252'); 
.libPaths(file.path(R.home(),'library')); 
library(doBy); 
options(timeout=600); 
library(plotrix); 
library(stringi); 
library(stringr); 
library(survival); 
library(rms); 
library(nnet); 
library(car); 
library(mgcv); 
pdfwd<-6; pdfht<-6; 
load('/Users/liukangning/Desktop/EmpowerROS/Analysis/cognitivedepression/datamergeall.Rdata'); 
if (length(which(ls()=='EmpowerStatsR'))==0) EmpowerStatsR<-get(ls()[1]); 
names(EmpowerStatsR)<-toupper(names(EmpowerStatsR)); 
originalVNAME<-names(EmpowerStatsR); 
ofname<-'cognitivedepression_15_tbl'; 
vname<-c(NA,'EDUCATION','EDUCATION.1','EDUCATION.2','EDUCATION.3','EDUCATION.4','EDUCATION.5','MARITAL','MARITAL.1','MARITAL.2','MARITAL.3','MARITAL.4','MARITAL.5','MARITAL.6','PIR','GENDER','GENDER.1','GENDER.2','AGE','RACE','RACE.1','RACE.2','RACE.3','RACE.4','RACE.6','RACE.7','BMI','BP','BP.1','BP.2','CORONARY','CORONARY.1','CORONARY.2','STROKE','STROKE.1','STROKE.2','DPQ','AFT','CERAD','DSST')[-1]; 
vlabel<-c(NA,'EDUCATION','  1','  2','  3','  4','  5','MARITAL','  Married','  Widowed','  Divorced','  Separated','  Never married','  Living with partner','PIR','GENDER','  male','  female','AGE','RACE','  Mexican American','  Other Hispanic','  Non-Hispanic White','  Non-Hispanic Black','  Non-Hispanic Asian','  Other Race - Including Multi-Racial','BMI','BP','  1','  2','CORONARY','  1','  2','STROKE','  1','  2','DPQ','AFT','CERAD','DSST')[-1]; 
varused4this <- c('EDUCATION','MARITAL','PIR','GENDER','AGE','RACE','BMI','BP','CORONARY','STROKE','DPQ','AFT','CERAD','DSST'); 
pkgs<-c('mgcv','geepack','lmtest'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(mgcv); 
library(geepack); 
library(lmtest); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'交互作用检验'; 
wd.subset=''; 
weights.var<- NA; 
yvname<-c('AFT','DSST','CERAD'); 
ydist<-c('gaussian','gaussian','gaussian'); 
ylink<-c('identity','identity','identity'); 
ylv<-c(0,0,0); 
colvname<-'STROKE'; colvlv<- 2; 
avname<-c('AGE','RACE','GENDER'); 
saf<-c(0,0,0); 
alv<-c(0,6,2); 
svname<-c('AGE','GENDER','RACE','MARITAL','PIR','EDUCATION','BMI','BP','CORONARY','STROKE'); 
sdf<-c(0,0,0,0,0,0,0,0,0,0); 
slv<-c(0,2,6,6,0,5,0,2,2,2); 
xvname<-c('DPQ'); 
xlv<-c(0); 
bvar<- NA; 
chk<- 1; 
cox<- 0; 
timevar<- NA; 
vname.start<- NA; 
subjvname<- NA; 
gee.TYPE<-NA; 
prn<-1; 
dec<-1; 
##R package## mgcv geepack lmtest ##R package##;
pvformat<-function(p,dec) {
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p)); colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  lw <- paste("<",substr("0.00000000000",1,dec+1),"1",sep="");
  pp[as.numeric(p)<(1/10^dec)]<-lw
  return(pp)
}
numfmt<-function(p,dec) {
  if (is.list(p)) p<-as.matrix(p)
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p));colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  pp[as.numeric(p)>10000000]<- "inf."
  pp[is.na(p) | gsub(" ","",p)==""]<- ""
  pp[p=="-Inf"]<-"-Inf"
  pp[p=="Inf"]<-"Inf"
  return(pp)
}
mat2htmltable<-function(mat) {
  t1<- apply(mat,1,function(z) paste(z,collapse="</td><td>"))
  t2<- paste("<tr><td>",t1,"</td></tr>")
  return(paste(t2,collapse=" "))
}
chkMatrixSinglar <- function(mat) {
  oo<-rep(0,ncol(mat))
  tmp<-mat[apply(is.na(mat),1,sum)==0,]
  oo[apply(tmp,2,function(z) length(levels(as.factor(z)))==1)]<-1
  rr<-cor(tmp)
  for (i in (2:ncol(rr))) if (sum(rr[1:(i-1),i]>0.999)>0) oo[i]<-1
  if (sum(oo)>0) print(rr)
  return(oo)
}
setgam<-function(fml,yi) {
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=gaussian(link="identity")))
  if (ydist[yi]=="binomial") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=binomial(link="logit")))
  if (ydist[yi]=="poisson") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=poisson(link="log")))
  if (ydist[yi]=="gamma") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=Gamma(link="inverse")))
  if (ydist[yi]=="negbin") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=negbin(c(1,10), link="log")))
  return(mdl)
}
setgee<-function(fml,yi) {
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") md<-try(geeglm(formula(fml),weights=wdtmp$weights,id=wdtmp[,subjvname],corstr=gee.TYPE,family="gaussian",data=wdtmp))
  if (ydist[yi]=="binomial") md<-try(geeglm(formula(fml),weights=wdtmp$weights,id=wdtmp[,subjvname],corstr=gee.TYPE,family="binomial",data=wdtmp))
  if (ydist[yi]=="poisson") md<-try(geeglm(formula(fml),weights=wdtmp$weights,id=wdtmp[,subjvname],corstr=gee.TYPE,family="poisson",data=wdtmp))
  if (ydist[yi]=="gamma") md<-try(geeglm(formula(fml),weights=wdtmp$weights,id=wdtmp[,subjvname],corstr=gee.TYPE,family="Gamma",data=wdtmp))
  if (ydist[yi]=="negbin") md<-try(geeglm.nb(formula(fml),weights=wdtmp$weights,id=wdtmp[,subjvname],corstr=gee.TYPE,data=wdtmp))
  return(md)
}
setgee2<-function(fmli,fmlo, yi, fml2) {
  id.gee=wdtmp[,subjvname]
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") {
    mdi<-try(geeglm(formula(fmli),id=id.gee,corstr=gee.TYPE,family="gaussian",weights=wdtmp$weights,data=wdtmp))
    mdo<-try(geeglm(formula(fmlo),id=id.gee,corstr=gee.TYPE,family="gaussian",weights=wdtmp$weights,data=wdtmp))
    if (!is.na(fml2)) md2<-try(geeglm(formula(fml2),id=id.gee,corstr=gee.TYPE,family="gaussian",weights=wdtmp$weights,data=wdtmp))
  }
  if (ydist[yi]=="binomial") {
    mdi<-try(geeglm(formula(fmli),id=id.gee,corstr=gee.TYPE,family="binomial",weights=wdtmp$weights,data=wdtmp))
    mdo<-try(geeglm(formula(fmlo),id=id.gee,corstr=gee.TYPE,family="binomial",weights=wdtmp$weights,data=wdtmp))
    if (!is.na(fml2)) md2<-try(geeglm(formula(fml2),id=id.gee,corstr=gee.TYPE,family="binomial",weights=wdtmp$weights,data=wdtmp))
  }
  if (ydist[yi]=="poisson") {
    mdi<-try(geeglm(formula(fmli),id=id.gee,corstr=gee.TYPE,family="poisson",weights=wdtmp$weights,data=wdtmp))
    mdo<-try(geeglm(formula(fmlo),id=id.gee,corstr=gee.TYPE,family="poisson",weights=wdtmp$weights,data=wdtmp))
    if (!is.na(fml2)) md2<-try(geeglm(formula(fml2),id=id.gee,corstr=gee.TYPE,family="poisson",weights=wdtmp$weights,data=wdtmp))
  }
  if (ydist[yi]=="gamma") {
    mdi<-try(geeglm(formula(fmli),id=id.gee,corstr=gee.TYPE,family="Gamma",weights=wdtmp$weights,data=wdtmp))
    mdo<-try(geeglm(formula(fmlo),id=id.gee,corstr=gee.TYPE,family="Gamma",weights=wdtmp$weights,data=wdtmp))
    if (!is.na(fml2)) md2<-try(geeglm(formula(fml2),id=id.gee,corstr=gee.TYPE,family="Gamma",weights=wdtmp$weights,data=wdtmp))
  }
  if (ydist[yi]=="negbin") {
    mdi<-try(geeglm.nb(formula(fmli),id=id.gee,corstr=gee.TYPE,weights=wdtmp$weights,data=wdtmp))
    mdo<-try(geeglm.nb(formula(fmlo),id=id.gee,corstr=gee.TYPE,weights=wdtmp$weights,data=wdtmp))
    if (!is.na(fml2)) md2<-try(geeglm.nb(formula(fml2),id=id.gee,corstr=gee.TYPE,weights=wdtmp$weights,data=wdtmp))
  }
  t.intr<-anova(mdi,mdo);
  print(fmli);print(fmlo);print(t.intr);
  p.intr<-pvformat(t.intr$"P(>|Chi|)",4)
  if (!is.na(fml2)) {
    print(fml2); t.intr2<-anova(md2,mdo); print(t.intr2);
    pintr<-paste(pintr," (",pvformat(t.intr2$"P(>|Chi|)",4)," #)",sep="");
  }
  return(list(mdi,mdo,p.intr))
}
setglm<-function(fml,yi) {
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") md<-try(glm(formula(fml),family="gaussian",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="binomial") md<-try(glm(formula(fml),family="binomial",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="poisson") md<-try(glm(formula(fml),family="poisson",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="gamma") md<-try(glm(formula(fml),family="Gamma",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="negbin") md<-try(glm.nb(formula(fml),weights=wdtmp$weights,data=wdtmp))
  return(md)
}
mdl2oo<-function(mdl, mdlo, xxname, opt, p.intr, mdlconti) {
  if (is.na(mdl[[1]][1])) return(list(rep("",times=length(xxname)),""))
  if (substr(mdl[[1]][1],1,5)=="Error") return(list(rep("",times=length(xxname)),""))
  if (is.na(mdlo[[1]][1])) return(list(rep("",times=length(xxname)),""))
  if (substr(mdlo[[1]][1],1,5)=="Error") return(list(rep("",times=length(xxname)),""))
  gs<-summary(mdl); print(mdl$formula); print(gs)
  gs2<-summary(mdlo); print(mdlo$formula); print(gs2)
  if (opt=="gam") {
    gsparm <- gs$p.table;tmpn<-gs$n; 
    pintr0<-anova.gam(mdl,mdlo,test="Chisq");print(pintr0)
    pintr<-pvformat(pintr0$"Pr(>Chi)"[2],4) 
    if (!is.na(mdlconti[[1]][1])) {
      pintr2<-anova.gam(mdlconti,mdlo,test="Chisq"); print(pintr2)
      pintr<-paste(pintr," (",pvformat(pintr2$"Pr(>Chi)"[2],4)," #)",sep="");
    }  
  } else {
    gsparm <- gs$coefficients;tmpn <- sum(gs$df[c(1,2)]);
    if (opt!="gee") {
      pintr0<-lrtest(mdl,mdlo); pintr<-pvformat(pintr0$"Pr(>Chisq)"[2],4);print(pintr0);
      if (!is.na(mdlconti[[1]][1])) {
        pintr2<-lrtest(mdlconti,mdlo); print(pintr2)
        pintr<-paste(pintr," (",pvformat(pintr2$"Pr(>Chisq)"[2],4)," #)",sep="");
      }  
    } else {pintr<-p.intr;}
  }
  gsp<-gsparm[match(xxname,rownames(gsparm)),]
  beta<-gsp[,1]; se<-gsp[,2]; pv<-gsp[,4]; ci1<- beta-1.96*se; ci2<- beta+1.96*se
  pvx<-substr(rep("****",length(pv)),1,(pv<=0.05)+(pv<=0.01)+(pv<=0.001)) 
  if (colprn==3) {pvv<-pvx;} else {pvv<-pvformat(pv,4);}
  if ((colprn!=2) & (gs$family[[2]]=="log" | gs$family[[2]]=="logit")) {
    o1<-paste(numfmt(exp(beta),dec)," (",numfmt(exp(ci1),dec),", ",numfmt(exp(ci2),dec),")",sep="")
  } else {
    if (colprn<3) {o1<-paste(numfmt(beta,dec), " (",numfmt(ci1,dec),", ",numfmt(ci2,dec),")",sep="")
    } else {o1<-paste(numfmt(beta,dec), "+",numfmt(se,dec),sep="");}
  }     
  o1<-paste(o1,pvv)
  return(list(o1,pintr))
}
mdl2bb<-function(mdl,xxname, opt) {
  gs<-summary(mdl); print(mdl$formula); print(gs)
  if (opt=="gam") {gsparm <- gs$p.table;} else {gsparm <- gs$coefficients;}
  bb<-round(gsparm[match(xxname,rownames(gsparm)),1],4)
  return(bb)
}
vv2fml<-function(vvname,vvlv,vsdf) {
  vvname_ <- vvname
  smoothvvi<-((vsdf=="s" | vsdf=="S") & vvlv==0)
  vvname_[vvlv>0]<-paste("factor(",vvname[vvlv>0],")",sep="")
  vvname_[smoothvvi]<-paste("s(",vvname[smoothvvi],")",sep="")
  fml<-paste(vvname_,collapse="+")
  return(fml)
}
addvia<-function(vvname,vvlv,vsdf,iavname) {    
  vvname_ <- vvname
  smoothvvi<-((vsdf=="s" | vsdf=="S") & vvlv==0)
  gamk<-(sum(smoothvvi)>0)
  vvname_[vvlv>0]<-paste("factor(",vvname[vvlv>0],")",sep="")
  vvname_[smoothvvi]<-paste("s(",vvname[smoothvvi],")",sep="")
  fml0<-paste(c(vvname_,paste("factor(",iavname,")",sep="")),collapse="+")
  vvname_2<-vvname
  vvname_2<-paste(vvname_2,"*factor(",iavname,")",sep="")
  vvname_2[smoothvvi]<-paste("s(",vvname[smoothvvi],", by=factor(", iavname, "))",sep="")
  vvname_2[vvlv>0]<-paste("factor(",vvname[vvlv>0],")*factor(",iavname,")",sep="")
  for (z in (1:length(vvname_))) {
    fml0<-c(fml0,paste(c(vvname_[-z],vvname_2[z]),collapse="+"))
  }
  vvia<-rep(FALSE,length(vvname))
  for (i in (1:ny)) {
    for (j in (1:nx)) {
        xxname<- xvianame_[[j]]; if (xlv[j]>0) xxname<-xxname[-1]
        xxname1<-paste(xxname,collapse="+")
        fml<-paste(yvname[i],"~",xxname1,"+",fml0); bb<-NA
        for (k in (1:length(fml))) {
          print(fml[k])
          if (gamk) {bb<-rbind(bb,mdl2bb(setgam(fml[k], i), xxname, "gam"));
          } else {
            if (!is.na(subjvname)) {bb<-rbind(bb,mdl2bb(setgee(fml[k], i), xxname, "gee"));
            } else {bb<-rbind(bb,mdl2bb(setglm(fml[k], i), xxname, "glm"));}
          }
        }
        bb<-bb[-1,]; tmprname_=paste("Model",1:length(fml)); rownames(bb)<-tmprname_; 
        tmpmdls_<-cbind(tmprname_,fml);  print(bb); print(tmpmdls_);
        bba<-(apply((bb< matrix(rep(bb[1,]-abs(bb[1,])*0.1,nrow(bb)),nrow=nrow(bb),byrow=TRUE)),1,sum)>0)
        bbu<-(apply((bb> matrix(rep(bb[1,]+abs(bb[1,])*0.1,nrow(bb)),nrow=nrow(bb),byrow=TRUE)),1,sum)>0)
        vvia <- (vvia | (bba | bbu)[-1])
    }  
  }
  #fmlx<-paste(c(paste("factor(",iavname,")",sep=""),vvname_[!vvia],vvname_2[vvia]),collapse="+")
  fmlx<-paste(c(vvname_[!vvia],vvname_2[vvia]),collapse="+")
  return(list(fmlx,vvia,vvname[vvia]))
}
colprn<-prn; if (colprn==0) colprn<-1;
if (!is.na(weights.var)) {weights<-WD[,weights.var];} else {weights<-1;}
WD<-cbind(WD,weights);
vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN]
sink(paste(ofname,".lst",sep=""))
Total_<-1; WD<-cbind(WD,Total_)
bvname<-bvar;
if (is.na(bvar)) {bvar<-"Total_"; bv<-WD[,bvar];bvb<-bvar;}
blv<-levels(factor(WD[,bvar])); nlb<-length(blv)
if (!is.na(bvname)) {
  bvb<-vlabelV[match(bvar,vnameV)][1];  if (is.na(bvb)) bvb<-bvar
  blb<-vlabelZ[match(paste(bvname,".",blv,sep=""),vnameZ)]
  blb[is.na(blb)]<-blv[is.na(blb)]
  stvlb<-c(paste(bvb,": ",blb,sep=""),"Total")
}
nx<-length(xvname); xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)]
ny<-length(yvname); yb<-vlabelV[match(yvname,vnameV)]; yb[is.na(yb)]<-yvname[is.na(yb)]
ns<-length(svname); sb<-vlabelV[match(svname,vnameV)]; sb[is.na(sb)]<-svname[is.na(sb)]
yref<-rep("Ref.",ny)
yref[ydist=="gaussian"]<-"0.0"
yref[ydist=="binomial"]<-"1.0"
yref[ydist=="poisson"]<-"1.0"
yref[ydist=="negbin"]<-"1.0"
clv<-levels(factor(WD[,colvname])); nlv<-length(clv)
colvb<-vlabelV[match(colvname,vnameV)]; if (is.na(colvb)) colvb<-colvname;
clb<-vlabelZ[match(paste(colvname,".",clv,sep=""),vnameZ)]
clb[is.na(clb)]<-clv[is.na(clb)]
xvname_ <- xvname
xvname_[xlv>0]<-paste("factor(",xvname[xlv>0],")",sep="")
xiaxname_ <- paste(xvname_,"+",xvname_,"*",colvname,sep="")
xvianame_<-list(NA); xvialbl_<-list(NA); 
for (j in (1:nx)) {
  tmpnewxv_<-NA; tmpnewxm_<-NA; tmpnewxb1<-""; tmpnewxb2<-"";
  if (xlv[j]==0) {
    for (k in (1:nlv)) {
      tmpnewxv_<-cbind(tmpnewxv_,WD[,xvname[j]]*(WD[,colvname]==clv[k]))      
      tmpnewxm_<-c(tmpnewxm_,paste(xvname[j],".",k,sep=""))
      tmpnewxb1<-c(tmpnewxb1,xb[j]);tmpnewxb2<-c(tmpnewxb2,clb[k])
    }
  } else {
    xjlv <- levels(factor(WD[,xvname[j]])); nxjlv<-length(xjlv)
    xjlb <- vlabelZ[match(paste(xvname[j],".",xjlv,sep=""),vnameZ)] 
    xjlb[is.na(xjlb)]<-xjlv[is.na(xjlb)]
    tmpxcc_<-rep(NA,nrow(WD)); tmpcode_<-0
    for (k in (1:nlv)) {
      for (z in (1:nxjlv)) {
        tmprzk_<-(WD[,xvname[j]]==xjlv[z] & WD[,colvname]==clv[k])
        tmpxcc_[tmprzk_]<-tmpcode_; tmpcode_<-tmpcode_+1
        tmpnewxv_<-cbind(tmpnewxv_,tmprzk_*1)
        tmpnewxm_<-c(tmpnewxm_,paste(xvname[j],".",z,".",colvname,".",k,sep=""))
        tmpnewxb1<-c(tmpnewxb1,xjlb[z]);tmpnewxb2<-c(tmpnewxb2,clb[k])
      }
    }
    tmpnewxv_<-cbind(tmpnewxv_,tmpxcc_);
    tmpnewxm_<-c(tmpnewxm_,paste(xvname[j],".",colvname,sep=""))
  }
  tmpnewxv_<-tmpnewxv_[,-1]; tmpnewxm_<-tmpnewxm_[-1]; 
  tmpnewxb1<-tmpnewxb1[-1]; tmpnewxb2<-tmpnewxb2[-1]
  colnames(tmpnewxv_)<-tmpnewxm_
  WD<-cbind(WD,tmpnewxv_)
  if (xlv[j]==0) {xvianame_[[j+1]]<-tmpnewxm_;} else {xvianame_[[j+1]]<-tmpnewxm_[-length(tmpnewxm_)];}
  xvialbl_[[j+1]]<-list(tmpnewxb1,tmpnewxb2)
}
xvianame_<-xvianame_[-1]; xvialbl_<-xvialbl_[-1]

allname<-c(yvname,xvname,colvname,bvar,avname,svname,subjvname)
allname<-allname[!is.na(allname) & allname>""]
if (!is.na(subjvname)) WD<-WD[order(WD[,subjvname]),]
wdtmp<-WD[(apply(is.na(WD[,allname]),1,sum)==0),]

fmlxx<-""; fmlzz<-""; fmltt<-"Crude"; fmlvv<-list(""); 
fmlgg<-ifelse(!is.na(subjvname), "gee", "glm");
if (length(avname)>0) {
  avb<-vlabelV[match(avname,vnameV)]; avb[is.na(avb)]<-avname[is.na(avb)]
  avb[((saf=="s" | saf=="S") & alv==0)]<-paste(avb[((saf=="s" | saf=="S") & alv==0)],"(smooth)")
  smoothav<-(sum((saf=="s" | saf=="S") & alv==0)>0)
  avbstr<-paste(avb,collapse=", ")
  fmlgg<-c(fmlgg,ifelse(!is.na(subjvname), "gee", ifelse(smoothav,"gam","glm")))
  fmlxx<-c(fmlxx,vv2fml(avname,alv,saf))
  fmlzz<-c(fmlzz,avbstr)
  fmltt<-c(fmltt,"Model I")
  fmlvv[[length(fmlvv)+1]]<-avname
  if (chk==1) {
    tmp<-addvia(avname,alv,saf,colvname)
    if (sum(tmp[[2]])>0) {
      fmlxx<-c(fmlxx,tmp[[1]]); 
      str1<- paste(vlabelV[match(tmp[[3]],vnameV)],collapse=", ")
      fmlzz<-c(fmlzz,paste(avbstr, " and the interaction terms for following variables:", str1));
      fmltt<-c(fmltt,"Model I*")
      fmlvv[[length(fmlvv)+1]]<-avname
      fmlgg<-c(fmlgg,ifelse(!is.na(subjvname), "gee", ifelse(smoothav,"gam","glm")))
    }
  }
}
if (length(svname)>0) {
  sb[((sdf=="s" | sdf=="S") & slv==0)]<-paste(sb[((sdf=="s" | sdf=="S") & slv==0)],"(smooth)")
  smoothsv<-(sum((sdf=="s" | sdf=="S") & slv==0)>0)
  svbstr<-paste(sb,collapse=", ")
  fmlgg<-c(fmlgg,ifelse(!is.na(subjvname), "gee", ifelse(smoothsv,"gam","glm")))
  fmlxx<-c(fmlxx,vv2fml(svname,slv,sdf))
  fmlzz<-c(fmlzz,svbstr)
  fmltt<-c(fmltt,"Model II")
  fmlvv[[length(fmlvv)+1]]<-svname
  if (chk==1) {
    tmp<-addvia(svname,slv,sdf,colvname)
    if (sum(tmp[[2]])>0) {
      fmlxx<-c(fmlxx,tmp[[1]]); 
      str1<- paste(vlabelV[match(tmp[[3]],vnameV)],collapse=", ")
      fmlzz<-c(fmlzz,paste(svbstr, " and the interaction terms for following variables:", str1));
      fmltt<-c(fmltt,"Model II*")
      fmlvv[[length(fmlvv)+1]]<-svname
      fmlgg<-c(fmlgg,ifelse(!is.na(subjvname), "gee", ifelse(smoothsv,"gam","glm")))
    }
  }
}
nmdl<-length(fmlxx)
if (nlb>1) {stmp_<-c(1:nlb,0);} else {stmp_<-0;}
oo1<-list(NA,NA,NA,NA,NA); oo2<-list(NA,NA,NA,NA,NA,NA,NA);
for (s in stmp_) {
  for (i in (1:ny)) {
    for (j in (1:nx)) {
      ifelse(xlv[j]>0, tmpvxz<-paste(xvname[j],".",colvname,sep=""), tmpvxz<-colvname) 
      wdtmp0<-WD[(apply(is.na(WD[,c(yvname[i],xvname[j],colvname)]),1,sum)==0),]
      tmpmean_<-tapply(wdtmp0[,yvname[i]],as.factor(wdtmp0[,tmpvxz]),mean)
      tmpnum_<-tapply(!is.na(wdtmp0[,yvname[i]]),as.factor(wdtmp0[,tmpvxz]),sum)
      if (ylv[i]>0) {
        tmpper_<-numfmt(tmpmean_*100,dec)
        stati<-paste(tmpnum_," (",tmpper_,")",sep="")        
      } else {
        tmpstd_<-tapply(wdtmp0[,yvname[i]],as.factor(wdtmp0[,tmpvxz]),sd)
        stati<-paste(numfmt(tmpmean_,dec),"+",numfmt(tmpstd_,dec),sep="")
      }
      numi<-tmpnum_;
      oo1<-rbind(oo1,list(s,i,j,numi,stati))
      for (k in (1:nmdl)) {
        tmpvn<-c(yvname[i],xvname[j],colvname);
        if (fmlvv[[k]][1]>"") tmpvn<-c(tmpvn,fmlvv[[k]])
        wdtmp<-WD[(apply(is.na(WD[,tmpvn]),1,sum)==0),]
        if (s>0) wdtmp<-wdtmp[wdtmp[,bvar]==blv[s],]
        nn<-tapply(!is.na(wdtmp[,yvname[i]]),as.factor(wdtmp[,tmpvxz]),sum)
        xxname<- xvianame_[[j]]; 
        if (xlv[j]>0) {
          xxname<-xxname[-1]
          fmli<-paste(yvname[i],"~",paste(xxname,collapse="+"),sep="")
        } else { 
          fmli<-paste(yvname[i],"~",paste(xxname,collapse="+"),"+factor(",colvname,")",sep="")
        }
        fmlo<-paste(yvname[i],"~",xvname_[j],"+factor(",colvname,")",sep="")
	    fml2<-paste(yvname[i],"~",xiaxname_[j],"+factor(",colvname,")",sep="")
        if (fmlxx[k]>"  ") {
          fmli<-paste(fmli,"+",fmlxx[k],sep="")
          fmlo<-paste(fmlo,"+",fmlxx[k],sep="")
	      fml2<-paste(fml2,"+",fmlxx[k],sep="")
        }
        if (s==0 & !is.na(bvname)) {
          fmli<<-paste(fmli,"+factor(",bvname,")",sep="")
          fmlo<<-paste(fmlo,"+factor(",bvname,")",sep="")
          fml2<<-paste(fml2,"+factor(",bvname,")",sep="")
        }
        mdlconti<-NA
        if (min(apply(wdtmp[,c(xxname,colvname)],2,function(m) length(unique(m))))>1) {
          if (fmlgg[k]=="gam") {
            if (nlv>2) mdlconti<-setgam(fml2,i)
            mdli<-mdl2oo(setgam(fmli,i),setgam(fmlo,i),xxname,"gam",NA,mdlconti)
          }
          if (fmlgg[k]=="gee") {
            if (nlv<=2) fml2<-NA
            tee<-setgee2(fmli,fmlo,i,fml2); mi<-tee[[1]]; mo<-tee[[2]]; p.intr<-tee[[3]]
            mdli<-mdl2oo(mi,mo,xxname,"gee",p.intr,NA)
          }
          if (fmlgg[k]=="glm") {
            if (nlv>2) mdlconti<-setglm(fml2,i)
            mdli<-mdl2oo(setglm(fmli,i),setglm(fmlo,i),xxname,"glm",NA,mdlconti)
          }  
        } else {
          mdli<-list(rep("",times=length(xxname)),"");
        }
        oo2<-rbind(oo2,list(s,i,j,k,mdli[[1]],mdli[[2]],nn))
      }
    }
  }
}

ncontx<-sum(xlv==0); ncatgx<-nx-ncontx
w<-""
if (ncontx>0) {
  tt1<-c("Model",paste(colvname,": ",clb,sep=""),"P interaction")
  nn1<-c("Y","Strata","X","Model",clb)
  for (i in (1:ny)) {
    if (ny>1) tt1<-rbind(tt1,c(paste("For ",yb[i],sep=""),rep(" ",nlv+1)))
    for (s in stmp_) {
       if (nlb>1) {
         ifelse(s>0, sslb<-stvlb[s], sslb<-stvlb[nlb+1])
         tt1<-rbind(tt1,c(sslb,rep(" ",nlv+1)))
       } else {sslb<-"Total";}
       for (j in (1:nx)) {
         if (xlv[j]==0) {
           if (nx>1) tt1<-rbind(tt1,c(xb[j],rep(" ",nlv+1)))
           for (k in (1:nmdl)) {
             oo2row <-(oo2[,1]==s & oo2[,2]==i & oo2[,3]==j & oo2[,4]==k)
             oo2row[is.na(oo2row)]<-FALSE
             mdloo<-oo2[oo2row,][[5]]; pia<-oo2[oo2row,][[6]];mdlnn<-oo2[oo2row,][[7]]
             tt1<-rbind(tt1,c(fmltt[k],mdloo,pia))
             nn1<-rbind(nn1,c(yb[i],sslb,xb[j],fmltt[k],mdlnn))
           }
         }
       }
    }
  }
}

if (ncatgx>0) {
  tt2<-c("Exposure",colvb,"N","Statistics")
  if (nmdl==1) mdlk<-1
  if (nmdl==2) mdlk<-2
  if (nmdl==3) {if (fmlxx[3]=="Model II") {mdlk<-c(2,3);} else {mdlk<-3;}}
  if (nmdl==4) {
    if (fmlxx[4]=="Model II") mdlk<-c(3,4)
    if (fmlxx[4]=="Model II*") mdlk<-c(2,4)
  }
  if (nmdl==5) mdlk<-c(3,5)
  mdlz<-(1:nmdl)[-mdlk]; nmdlz<-length(mdlz)
  tt2<-c(tt2,fmltt[mdlk])
  if (nmdlz>0) {tt3<-c("Exposure",colvb,fmltt[mdlz]);} else tt3<-NA
  tbncol<-length(tt2)
  nn2<-c("Y","Strata","X","Model","Nobs","Tot N")
  for (i in (1:ny)) {
    if (ny>1) {
      tt2<-rbind(tt2,c(paste("For ",yb[i],sep=""),rep(" ",tbncol-1)))
      if (nmdlz>0) tt3<-rbind(tt3,c(paste("For ",yb[i],sep=""),rep(" ",nmdlz+1)))
    }
    for (s in stmp_) {
       if (nlb>1) {
         ifelse(s>0, sslb<-stvlb[s], sslb<-stvlb[nlb+1])
         tt2<-rbind(tt2,c(sslb,rep(" ",tbncol-1)))
         if (nmdlz>0) tt3<-rbind(tt3,c(sslb,rep(" ",nmdlz+1)))
       } else {sslb<-"Total";}
       for (j in (1:nx)) {
         if (xlv[j]>0) {
           oo1row <-(oo1[,1]==s & oo1[,2]==i & oo1[,3]==j)
           oo1row[is.na(oo1row)]<-FALSE
           mdlcc<-cbind(c(xb[j],xvialbl_[j][[1]][[1]]),c(" ",xvialbl_[j][[1]][[2]]))
           if (nmdlz>0) mdlcx<-mdlcc
           mdlcc<-cbind(mdlcc,c(" ",oo1[oo1row,][[4]]),c(" ",oo1[oo1row,][[5]]))
           pia<-c("P interaction"," "," "," ")
           if (nmdlz>0) pix<-c("P interaction"," ")
           for (k in mdlk) {
             oo2row <-(oo2[,1]==s & oo2[,2]==i & oo2[,3]==j & oo2[,4]==k)
             oo2row[is.na(oo2row)]<-FALSE
             mdloo<-c("","Ref.",oo2[oo2row,][[5]]); 
             pia<-c(pia,oo2[oo2row,][[6]]);
             mdlnn<-oo2[oo2row,][[7]]; sumnn<-sum(mdlnn)
             mdlcc<-cbind(mdlcc,mdloo)             
             nn2<-rbind(nn2,c(yb[i],sslb,xb[j],fmltt[k],paste(mdlnn,collapse=";"),sumnn))
           }
           tt2<-rbind(tt2,mdlcc,pia)
           if (nmdlz>0) {
             for (k in mdlz) {
               oo2row <-(oo2[,1]==s & oo2[,2]==i & oo2[,3]==j & oo2[,4]==k)
               oo2row[is.na(oo2row)]<-FALSE
               mdloo<-c("","Ref.",oo2[oo2row,][[5]]); 
               pix<-c(pix,oo2[oo2row,][[6]]);
               mdlnn<-oo2[oo2row,][[7]]; sumnn<-sum(mdlnn)
               mdlcx<-cbind(mdlcx,mdloo)             
               nn2<-rbind(nn2,c(yb[i],sslb,xb[j],fmltt[k],paste(mdlnn,collapse=";"),sumnn))
             }
             tt3<-rbind(tt3,mdlcx,pix)
           }
         }
       }
    }
  }
}
sink()
prnopt<-c("β (95%CI) Pvalue / OR (95%CI) Pvalue", "β (95%CI) Pvalue", "β+se / OR (95%CI) *P<0.05 **P<0.01 ***P<0.001")
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,paste("<h2>", title, "</h2>"))
 

if (ncontx>0) {
  w<-c(w,"</br><table border=3>", mat2htmltable(tt1), "</table>")
  w<-c(w,"</br>Results in table:", prnopt[colprn])
  if (nlv>2) w<-c(w,"</br># interaction trend test (treat ", colvb, " as continuous)")
  w<-c(w,paste(c("</br>结果变量:",yb),collapse=" "))
  w<-c(w,paste(c("</br>危险因素:",xb[xlv==0]),collapse=" "))
  w<-c(w,paste(c("</br>效应修饰因子:",colvb),collapse=" "))
  for (k in (1:nmdl)) {
    if (fmlxx[k]>"") w<-c(w,paste("</br>",fmltt[k]," 调整了: ",fmlzz[k],sep=""))
  }
  if (sum(fmlgg=="gee")>0) w<-c(w, paste("</br>Generalized estimate equation were used, subject ID=", subjvname, "(", gee.TYPE,")",sep=""))
  w<-c(w,"</br></br>各模型所用的样本量</br><table border=3>", mat2htmltable(nn1), "</table>")
}
if (ncatgx>0) {
  w<-c(w,"</br><table border=3>", mat2htmltable(tt2), "</table>")
  w<-c(w,"</br>Results in table:", prnopt[colprn])
  if (nlv>2) w<-c(w,"</br># interaction trend test (treat ", colvb, " as continuous)")
  w<-c(w,paste(c("</br>结果变量:",yb),collapse=" "))
  w<-c(w,paste(c("</br>危险因素:",xb[xlv>0]),collapse=" "))
  w<-c(w,paste(c("</br>效应修饰因子:",colvb),collapse=" "))
  for (k in mdlk) {
    if (fmlxx[k]>"") w<-c(w,paste("</br>",fmltt[k]," 调整了: ",fmlzz[k],sep=""))
  }
  if (!is.na(tt3[1])) {
    w<-c(w,"</br><table border=3>", mat2htmltable(tt3), "</table>")
    w<-c(w,"</br>Results in table:", prnopt[colprn])
    w<-c(w,paste(c("</br>结果变量:",yb),collapse=" "))
    w<-c(w,paste(c("</br>危险因素:",xb[xlv>0]),collapse=" "))
    w<-c(w,paste(c("</br>效应修饰因子:",colvb),collapse=" "))
    for (k in mdlz) {
      if (fmlxx[k]>"") w<-c(w,paste("</br>",fmltt[k]," 调整了: ",fmlzz[k],sep=""))
    }
  }
  if (sum(fmlgg=="gee")>0) w<-c(w, paste("</br>Generalized estimate equation were used, subject ID=", subjvname, "(", gee.TYPE,")",sep=""))
  w<-c(w,"</br></br>各模型所用的样本量</br><table border=3>", mat2htmltable(nn2), "</table>")
}
w<-c(w,wd.subset)
w<-c(w,paste("</br>此表用易侕统计软件 (www.empowerstats.com) 和R软件生成，生成日期：",Sys.Date()))
w<-c(w,"</body></html>")
fileConn<-file(paste(ofname,".htm",sep="")); writeLines(w, fileConn)